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SUMMARY 


The governing differential equations of motion for a high speed 
cylindrical roller bearing are programmed for numerical solution and 
plotted output. This computer program has the capability of performing 
a two-dimensional or three-dimensional simulation. Two numerical solu- 
tions of the governing differential equations were obtained to simulate 
the motion of a roller bearing, one for the two-dimensional equations of 
motion and one for the three-dimensional equations of motion. Computer 
generated plots were obtained and present such data as roller/cage 
interaction forces, roller /race traction forces, roller /race relative 
slip velocities and cage angular speed over a non-dimensional time 
equivalent to 1.2 revolutions of the inner race. Roller axial displacement, 
roller skew angle and skew moment are also plotted for the three-dimensional 
solution. The trajectory of the cage center is plotted for both the two- 


and three-dimensional solutions. 



INTRODUCTION 


The prediction of transient dynamic behavior of high speed cylindrical 
roller bearings is necessary to evaluate new design concepts and design 
changes in existing bearings. Many phenomena, such as cage-roller impacts 
and roller-flange impacts, which limit bearing life or affect bearing per- 
formance are of a transient nature and may not be predicted using steady 
state simulations. The prediction of roller or cage stability must be 
done using methods of transient dynamic analysis with the proper nonlinear 
models. Part I [1] of this report presents the derivation of the equations 
of motion necessary to do the transient dynamic analysis of roller bearings. 

Over the last decade much attention has been focused on the problem 
of rolling element bearing transient analysis. Walters [2] was one of the 
first investigators to attempt the solution of the total bearing dynamics 
problem with realistic lubrication models. He made the basic assumption 
in his analysis that the ball-race elasticity analysis can be performed 
independently prior to the dynamics analysis because the ball-race normal 
forces are large and the other forces, i.e. traction, ball spinning, have 
little influence on the ball-race force balance. This same assumption was 
used in [1] and is programmed in this report. Walters used a fourth order 
Runge-Kutta scheme to integrate the equations of motion in contrast to the 
use of the multivalue method of Gear [3] for stiff systems which is used in 
this report. 

Gupta [4] has presented a series of papers which treat the dynamics of 
both ball and cylindrical roller bearings. The force models utilized by Gupta 
are essentially the same as used in this computer program. The formula- 
tion of the differential equations, as programmed in this report, is different 
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from [4] and is discussed in [1], A Runge-Kut ta-Mer son method was used 
to numerically integrate the differential equations of motion in [4] again 
in contrast to the use of the multivalue method of Gear [3] in this report. 

The purpose of this research program was to develop computationally 
efficient models which describe the interactions among forces, relative 
velocities and relative displacements in a high speed lightly loaded cylin- 
drical roller bearing and to program the resulting differential equations 
of motion for a digital computer to obtain numerical solutions to these dif- 
ferential equations. These objectives have been achieved and the resulting 
computer program is presented in subsequent sections of this report. 

The section ORGANIZATION OF THE COMPUTER PROGRAM discusses each subroutine. 
The dependent variables are defined together with their ordering scheme. 

The results of a two-dimensional and a three-dimensional simulation 
are presented in RESULTS to illustrate the use of the computer program. 

These output data are in the form of plots of quantities of interest such 
as roller-race traction forces and roller-cage impact forces. All the 
plotted output data are discussed in DISCUSSION OF RESULTS. The final 
section contains the CONCLUDING REMARKS. The Appendices contain detailed 
information on the input data and output data. 
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ORGANIZATION OF THE COMPUTER PROGRAM 


The purpose of the computer program is to obtain a numerical solution 
to the governing differential equations of motion described in Part I [1] . 
These equations are represented as a system of first order nonlinear dif- 
ferential equations of the form 

= 1 (Z. T) (1) 

where ^ is the vector of state variables 
T is the dimensionless time 
^ is the vector function 

The computer program is organized into three major parts: the m.ain 

program which controls the input and output, the differential equation 
solver, and the subroutines which are used to evaluate the vector function 
f T) . Fig. 1 shows the interrelationship of all the subroutines in the 

computer program. 

The multivalue method of Gear [3] was utilized in the development of 
the computer program. This method has been shown to be effective in the 
solution of stiff systems of equations. The predictor step is taken by ex- 
trapolating out to a new value of the state vector. The extrapolation func 
tion is a weighted combination of values of the state vector at previous 
times. The corrector step is based on the implicit method of solving dif- 
ferential equations. 

The multivalue method proved to be successful and necessary when the 
problem of the roller impacting the cage and flange was introduced. The 
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abrupt change in compliance upon impact causes the eigenvalues of the 
differential equations to change abruptly. These eigenvalues can range 
over several orders of magnitude when some rollers are under impact and 
others are not. 

Each major subroutine in the computer program will now be discussed. 
Program ILLIBRG 

The purpose of the main program is to control the input and the out- 
put of the program and to call the differential equation solver. This pro- 
gram also calculates constants for use throughout the program and converts 
input quantities to other units which are more convenient for use in the 
program. The maximum and minimum allowable step sizes together with the 
specification of the total time over which the equations are integrated are 
controlled from this program. 

The input is read in card-image form and includes the geometrical and 
inertial properties of the roller and the cage, the properties of the lubri- 
cant and the initial conditions of the bearing state vector y. The initial 
conditions on the first derivatives with respect to time are normally read in 
units of the epicyclic speed and are converted in the computer program to 
units of the speed as given in equation (73) of Part I [1] . 

The output of the program is the state vector which is printed at 
time intervals subject to the user*s option. Other data, including the rol- 
ler traction forces, the relative slip velocities, the roller skewing moments 
and the trajectory of the cage center, are plotted from subroutine PLODAT, 
which is the most convenient way to present the large amount of output data. 

In addition to the plotting subroutine, files have been set up to store the 
output data of interest at a time interval set by the user. After the program 
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has terminated, the data stored on these files can be used as data in a 
special purpose plotting program to machine plot the data, or as input 
data for a subsequent run. 

Differential Equation Solver 

The differential equation solver used was the multivalue value 
method for stiff systems programmed in GEAR. This subroutine is called 
by the main program with a desired step size. The subroutine GEAR and its 
subsidiary subroutines perform the integration for the desired step size 
and returns to the main program where the time variable is updated. The 
subroutine GEAR used in this program was developed at the Lawrence Livermore 
Laboratory [5] . 

Subroutine for Plotting Data 

Operational data generated by the computer program is plotted according 
to the instructions in subroutine PLODAT. The data plotted by this sub- 
routine are the net circumferential forces on the roller, the traction forces 
at the inner and outer races, the roller to race relative slip velocities 
at the inner and outer race, the ratio of the cage angular velocity to the 
epicyclic speed and the trajectory of the cage center. Additionally, for 
the three dimensional problem, the roller axial displacement, the roller 
skew and the skew moment due to traction force distribution between the 
roller and the races are plotted. 

Title blocks are inserted on each plot and are specific to the problem 
being solved. The plot statements used in PLODAT are specific to the Illi- 
nois Calcomp'^ Graphics system and would have to be modified to 
fit the plotting systems of other computer installations. 
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The plot data are stored in arrays in the main program, ILLIBRG, 
in two dimensional array form. After the integration has been completed, 
these data are rearranged and packed into one dimensional arrays with two 
blank spaces separating the data for each roller. These spaces are used 
to store the scaling factors needed to plot the data in subroutine PLODAT. 

Subroutine for the Differential Equation Solver 

Subroutine GEAR requires a subroutine to calculate the vector function 
f^ (^, T) in equation (1). This subroutine is named DIFFUN. 

Subroutine DIFFUN calculates many of the forces and torques acting in 
the roller bearing system. All the remaining subroutines are called di- 
rectly or indirectly by DIFFUN in the calculation of the vector function 
1 (X, T), Fig. 1. 

Subroutine SHVEL 

Subroutine SHVEL takes the user supplied values of total shaft linear 
deflection, the shaft rotation speed, the total angle of misalignment be- 
tween the inner race and the outer race and whether the linear deflection 
and/or the angular misalignment is rotating or steady. This subroutine 
returns values of the linear deflections of the shaft in the x and y direc- 
tions and the angular misalignments measured about the x and y axes as a 
function of time. The overall computer program, as it is now written, can 
only utilize SHVEL if the shaft speed is constant. If studies are to be 
made for an accelerating bearing, the equation 

®shaft v/ “shaft (2) 
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would have to be solved using GEAR and the value of 0 , would be 

shaft 

required as input to SHVEL together with the time varying value of 
Function WSF 

The function WSF is designed to give the shaft angular speed as a 
function of time. Presently the function is programmed to have a con- 
stant speed with time, but any speed-time relationship may be programmed. 

If the effect of a variation of shaft speed is desired, equation (2) 
would have to be programmed in DIFFUN and the result, would have 

to be included as an input parameter to SHVEL. 

Subroutine RQLDIS 

Subroutine ROLDIS is called by ILLIBRG and DIFFUN. The method of 
Harris [6] and Liu [7] is programmed in this subroutine and its subsidiary 
routines. The geometry of the roller and the external displacements of 
the inner race with respect to the outer race are inputs to this sub- 
routine. ROLDIS then organizes this information and returns with the static 
load distribution on each of the rollers. 

Subroutine SETROL 

SETROL computes the deviations of the roller surface from a pure cyl- 
inder due to crowning. The vector of variables, C^, (refer to equation 
(8), [1]) at the center of each slice is an output of this subroutine. 

The roller centrifugal force per unit cage speed (in rpm) squared is also 
calculated in this subroutine. 

Subroutine SETVEC 

SETVEC calculates the values of the sine, cosine and the sign of the 
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cosine of the angles corresponding to the positions of the rollers in the 
bearing. 

Subroutine INITL 

This subroutine provides an initial set of values for the and 

the 6°^'s in the solution of the nonlinear algebraic system of equations 
in ROLDIS. This subroutine is only called once; on the second and later 
calls, the previous solutions are used as the initial values for the iter- 
ation scheme. 

Subroutine SETMAT 

SETMAT sets up the Jacobian matrix for the Newton-Raphson iteration 
for the roller load distribution calculation. This routine is called at 
every iteration. 

Subroutine NRHAP 

This subroutine controls the Newton-Raphson iteration for the load 
distribution on the rollers. The systems of (2 X 2) equations are solved 
and the resulting incremental changes of the dependent variables are com- 
pared with the lower limit. If the changes are less than 10 the itera- 
tion process is assumed to have converged. 

Function COEFF 

The correlation equation which relates the rolling velocity, sliding 
velocity, contact pressure and viscosity of the lubricant to the parameters 
necessary to obtain the traction coefficient is programmed in this function 
subroutine. The traction coefficient for the oil film between the roller 
and the races is the output of COEFF. The equations (20) - (23) of Part I 
[1] are programmed in COEFF. 
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Subroutine RLCGF 


The regression equations for the dimensionless normal forces be- 
tween the roller and cage are programmed in RLCGF. These equations are 
taken from Dowson, Markho and Jones [8] • 

Subroutine FRICTF 

The regression equations from [8] for the dimensionless tractive 
forces between the roller and the cage are programmed in FRICTF. The 
two outputs of this subroutine are the dimensionless tractive force on 
the roller and the dimensionless tractive force on the cage. 

Subroutine R1245 

The drag on the cylindrical surface of the roller is calculated in 
R1245. The output of this subroutine is the drag torque on the outer 
surface of the cylinder, The expressions used in this subroutine 

are based on the models used by Rumbarger, et. al. [10] . 

Subroutine ROLEND 

The drag torques on the ends of the roller as the roller turns about 
its own axis are calculated in ROLEND. The expressions in this subroutine 
are based on the models used by Rumbarger, et. al. [10]. 

Subroutine CGTRQ 

The drag torque on the outer surface and sides of the cage are cal- 
culated in CGTRQ. The area of the outer surface of the cage is estimated 
to be the surface area of the equivalent cylinder less the area of the 
pockets for the rollers. The expressions in this subroutine are based 
on the models in [10] . 
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Subroutine LANDLIP 


This subroutine calculates the horizontal driving force and the 
driving torque on the roller as a result of shearing the lubricant be- 
tween the roller end and the inner race guide flange. The axial clear- 
ance between a roller end and the guide flange is assumed to be equal on 
both sides for purposes of this calculation for both the two- and three- 
dimensional simulations. The outputs of this subroutine are a net driving 
torque on the roller Tj^lMD ^ driving force F^land* 

Subroutine GQU3Z 

This subroutine performs a composite Legendre-Gauss quadrature using 
a 96-point formula. The function which is numerically integrated is the 
pressure between the cage and the inner race guide flanges due to hydro- 
dynamic action. The output of this subroutine is net force on the cage 
due to this pressure. 

Subroutine CGPRES 

The pressure function which is utilized in GQU3Z is evaluated in 
this subroutine. A check is made in this subroutine to assure that only 
the positive pressures are integrated. 

The Dependent Variables 

The dependent variables in the computer program are in the array 
Y(J). The number of dependent variables for the two dimensional problem 
is 3(NZ)+6 where NZ is the number of rollers in the bearing. The first 
six elements of Y(J) are used to describe the state variables of the cage. 
The first two elements are the dimensionless velocity and displacement of 
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the cage in the x direction. The next two elements are the dimensionless 
velocity and displacement of the cage in the y direction. The fifth and 
sixth elements are the dimensionless angular velocity and angular dis- 
placement of the cage. 

The next group of three variables (seventh, eighth and ninth) repre- 
sents the dimensionless relative velocity and displacement between the 
cage and the first roller and the angular velocity of the first roller 
about its own axis. The following group of three variables is associated 
with the second roller and so forth, until the state variables of all NZ 
rollers have been described. 

The three dimensional problem requires an additional 3(NZ) variables 
which are concatenated to the original 3(NZ)+6 variables described above. 
The first group of three variables represents the skewing angular velocity, 
the angular displacement and the axial displacement of the first roller. 

The next (NZ -1) groups of three variables represent corresponding quan- 
tities for the remaining (NZ -1) rollers. 



RESULTS 


The roller bearing dynamics program was utilized to illustrate the 
dynamic behavior of an inner race guided high speed lightly loaded cylin- 
drical roller bearing. The geometry of the bearing used in this investi- 
gation is given in Table 1. The material of the rollers, races and the 
cage is a conventional bearing steel. The relevant material properties 
of the steel and the lubricant together with selected inertial properties 
of the cage are given in Table 2, A steady state condition was obtained for 
a two-dimensional simulation at a shaft angular speed of 70,000 revolutions 
per minute with no misalignment and a rotating linear deflection of the shaft 
of 0.020 mm, the state variables of which were used as the initial conditions 
for the two cases analyzed in this report. These initial conditions on the 
first derivatives with respect to time were non-dimens ionali zed with 

respect to the theoretical epicyclic speed eqn, (3), 


epx shaft 

2 (1 H- f) 


(3) 


is the theoretical epicyclic speed (radian/s) 
is the shaft angular velocity (radian/s) 

is the roller radius (mm) 

is the radius of the inner race (ram) 

The two-dimensional simulation was run for 1,2 revolutions of the inner 
ring. The loading condition was a rotating linear deflection of the shaft 
of 0.025 mm with no misalignment. The three-dimensional simulation was also 
run for 1,2 revolutions of the inner ring with a loading condition of a 

bearing misalignment of 0.002 radians and a steady shaft deflection in the y 

direction of 0,020 mm. 


where co 

epx 

shaft 

r 

R 
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The data obtained for the above two conditions are presented in Figures 

2 through 18, The data presented for the two-dimensional case are the trac- 
tion forces at the inner race and the outer race, the roller slip velocities 
at the inner race and outer race, the net circumferential forces between the 
rollers and the cage, the cage angular speed normalized on the theoretical 
epicyclic speed, and the trajectory of the cage center normalized on the 
radial clearance. The lines in the figures with multiple plots are labelled 
A through G corresponding to the seven rollers in the bearing. 

The data presented for the three-dimensional case include the traction 
forces at the inner race and outer race, the net circumferential forces be- 
tween the roller and the cage, the roller slip velocities at the inner race 
and outer race, the cage angular speed normalized on the theoretical epicyclic 
speed, the roller skew angle, the roller skew moment due to misalignments 
of the roller, and the dimensionless roller axial displacement. The input 
data corresponding to the two conditions described above are given in Tables 

3 and 4. Appendix A, INPUT DATA, gives a detailed explanation of all the 


elements in these data sets 


TABLE 1 


BEARING GEOMETRY 


(ram) 

Bore 25 

Inner Race Diameter 31 

Guide Flange Height 2 

Corner Radius, Inner Race Guide Flange 0.762 

Outside Diameter 52 

Outer Race Diameter 45 

Roller Diameter 7 

Roller Length 7 

Crown Radius 635 

Corner Radius of Roller 0,762 

Cage Outside Diameter 42 

Cage Inside Diameter 36 

Cage Width 13 

Roller/Cage Pocket Clearance, Circumferential 0,406 

Roller/Cage Pocket Clearance, Axial 0.229 

Pitch Diameter (Reference) 38 

Diametral Clearance 0,025 

Cage/Outer Race Diametral Clearance 4 

Cage/Inner Race Guide Flange Diametral Clearance 0.356 

Width of Guiding Lands/Cage Surfaces (2) 3 

Roller /Inner Race Guide Flange Clearance 0,025 
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TABLE 2 


MATERIAL, INERTIAL AND OPERATIONAL PARAMETERS 


Bearing 

Modulus of Elasticity; Roller, Gage, Races 
Poisson’s Ratio; Roller, Cage, Races 
Roller Weight 

Density of Roller Material 

Polar Moment of Inertia of the Cage 

Cage Weight 

Ratio of Oil Volume/Total Volume in Bearing 
Number of Rollers 


2.0 X lO^^N/m^ 

0.30 

0.020 N 

7.75 X 10^ kg/m^ 
7.71 X 10"^ kg-m^ 
0.176 N 
0.20 
7 


Lubricant (Mil-L-7808) 

Absolute Viscosity 
Density 


4.152 cp 
830 kg/m^ 
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TABLE 3 


IMPaT DATA »HD 
PRESENTED IN 


INITIAL CONDITIONS FOR 2-*D SIHDLATION 
THE FORMAT DESCRIBED IN APPENDIX A 


7 8 2 

30000000. 30000000. 
1.232 .2756 

.2756 25.0 

.68213900000000B-04 
4.151985 .03 

0.00 0.00 
70000.00 
.0002 0.00 
-01 

-.9874705645D-0 1 
-*.2213403995D+00 
2459280210D+00 
. 1207631242D+01 
.3836438428D+00 
.3601629631D-01 
-.4948767425D+00 
.1092699572D+00 
. 10757 14868D+00 
.2 157924327D+03 


30000000. .30 .30 .30 

.001 .0005 .030 .030 

.004603 .72538000000000E-03 
1.398 .03961 .11470 .007 

.20 

.001 1-0 

1.20 .05 


.8191528305D+00 
.8637317610D+00 
.370477201 OD+00 
. 49195H1581D+00 
--1730080141D+00 
-.2337202076D+00 
-4623993970D+00 
-.1780332250D+00 
.4433184912D+00 
.4899629083D+04 


.4095124256D400 

-3615456008D+01 

.4696269451D+01 

.48650912SOC401 

.4810836815D+01 

.4908954248D+01 

.4801317850D+01 

.4737511335D401 

.4737316664E401 


.016 


.12 


TABLE 4 


INPUT DATR AND 
PRESENTED IN 


INITIAL CONDITIONS FOR 3-D SIMULATION 
THE FORMAT DESCRIBED IN APPENDIX A 


7 8 3 

30000000. 30000000. 
1.232 .2756 

.2756 25.0 

.68213900000000E-04 
4.151985 .03 

0.00 0.00 
70000 .00 
-0002 0.00 
-01 

-. 987470564 7D-01 
-.2213403995D+00 
-. 2459280210D+00 
. 1207631242D401 
.3836438428D+00 
. 3601629631D-01 
--4948767426D400 
.1 092699572D+00 
. 10757 148680*00 
0 . 

0 - 

0 . 

0 . 

0 . 

0 . 

0 - 

. 1328643991D403 


30000000. .30 .30 -30 

.001 -0005 .030 -030 

.004603 -72538000000000E-03 
1.398 .03961 .11470 .007 

.20 

-001 1.0 


1-20 .05 

.8191528305D400 
.86373 1761 1D+00 
.3704772010D+00 
.49195415810*00 
-.17300801410*00 
--2337202076D400 
.46239939700*00 
-.17803322500*00 
.44331849120+00 
0 . 

0 - 

0 . 

0 . 

0 . 

0 . 

0 . 

.3844584514D+04 


.40951242571*00 

0 . 

.4696269452D+01 
.4865091250D+01 
.4810836815D+01 
-4908954249t*01 
. 4801317851E+01 
.47375 11335D+01 
.47373166650*01 
0 . 

0 . 

0 . 

0 . 

0 . 

0 . 

0 . 


.016 


.12 
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DISCUSSION OF RESULTS 


The results presented in this report are intended to illustrate the 
behavior of the computer program and to show the type of output data which 
may be obtained. The bearing used in these illustrative studies represents 
a high speed cylindrical roller bearing with the exception that the number 
of rollers was reduced from fourteen to seven. This was done to conserve 
computer time, as the primary purpose of this project was to develop a working 
computer program. 

The effects of heat transfer were not considered in this project nor 
were they programmed in the computer code. The lubricant properties are 
specified as input data and they correspond to an operating temperature of 
the bearing. If a lubricant other than Mil-L-7808 oil is to be used, the 
function subprogram COEFF will have to be changed to provide the traction 
coefficient for the new lubricant, 

Two-Dimensional Simulation 

The two-dimensional simulation initial conditions and loading param- 
eters represent the effect of a suddenly applied unbalance condition in 
an already unbalanced shaft. The maximum roller load increased from 591 N 

to 960 N. The effect of this sudden loading is shown in Figures 2 through 8. 

The dimensionless normal cage forces acting on the rollers, FC, are 
plotted in Figure 2, The purpose of this plot is to show the time and lo- 
cation of cage roller impacts, A positive value of FC indicates that the 
cage is pushing the roller while a negative value indicates that the roller 
is pushing the cage. For most of the 1,2 inner ring revolutions the value of 

FC appears to be zero, but in fact is not. The components of FC are the normal 

approach hydrodynamic models and the Hertzian impact models. The magnitudes 
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Figure 2. The Dimensionless Normal Cage Forces Acting on the Rollers for the 2-D Simulation Versus Time 
Measured in Revolutions of the Inner Race at 70,000 RPM 



of the hydrodynamic forces are very small relative to those forces gen- 
erated during Hertzian impact. When the data are scaled for plotting, 
the hydrodynamic forces will be plotted close to zero. 

To examine the magnitude of the impact forces, consider the force 
between roller C and the cage at about 1,07 revolutions. The maximum 
impact force is (,191) (960) = 183 N (41 Ib^), This force corresponds 
to a maximum Hertz stress between the cage pocket and roller C of 5,19 x 
10^ Pa (75,225 psi) , 

The dimensionless traction forces on the rollers at the inner race, 
FIN, are plotted in Figure 3. The traction force is directly propor- 
tional to the normal force at the contact, which accounts for the periodic 
shapes as shown. In a high speed bearing, the rollers are usually pressed 
against the outer race due to the action of centrifugal force. Depending 
on the amount of diametral clearance, there is usually no contact between 
the roller and the inner race except in the load zone of the bearing. The 
time required for a roller to traverse this load zone is 0,44 revolutions 
of the inner race. 

The magnitudes of the inner race/ roller traction forces are small, 
about 7 N (1,58 Ib^) , due to the low normal forces between the rollers 
and the inner race. The dimensionless outer race/roller traction forces, 
Four, however, achieve larger magnitudes. Figure 4, The magnitudes of 
these forces reach 48 N (10.8 Ib^) . 
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70000RPM SHRFT DEFLECTI 3N=0. 025MM MISRLIGNMENT=O.OOOORflO. 2-D MOTION 
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Figure 4a. The Dimensionless EHD Traction Forces on the Rollers at the Outer Race for the 2~D Simulation 
Versus Time Measured in Revolutions of the Inner Race at 70,000 RPM 



K) 

-O' 


70000RPM SHAFT DEFLECT 1 0N = O. 025MM MISRL IGNMENT = 0. OOOORflD. 2-D MOTI0N 
ROTATING DEFLECTION STEADY MISALIGNMENT MAX. ROLLER LORD IN) = 960. 



Figure 4b. The Dimensionless EHD Traction Force on Roller B at the Outer Race for the 2-D Simulation 
Versus Time Measured in Revolutions of the Inner Race at 70,000 RPM 




In the high speed bearing, the centrifugal force is significant and 
accounts for a large magnitude normal force between the rollers and the 
outer race around the circumference of the bearing. Since traction force 
is proportional to the normal force, the large magnitude of the traction 
force should be expected. Note that many of these dimensionless traction 
forces. Four, hover near zero (0 to 0.001). Note that the rollers under- 
going the largest changes in traction forces are those rollers going through 
the load zone. This correlates with the results in Figure 3, 

The relative slip velocities between the roller and the inner race, 
USIN, are plotted in Figure 5, Note that all the rollers have slip veloci- 
ties between 12 and 16 m/s. These appear to be approaching an equilibrium 
condition and correspond to the traction forces required to equilibrate 
the rollers because of the small inner race/roller normal forces. 

The relative slip velocities, USOUT, between the rollers and the 
outer race are shown in Figure 6, Note that all the slip velocities tend 
toward zero. This correlates with the values of the traction forces in 
Figure 4. The smallest deviation of relative slip velocity from zero is 
sufficient to cause a significant change in the traction forces FOUT, The 
non-zero values of USOUT in the first 0.10 revolutions of the inner race 
correspond to the set of initial conditions. At equilibrium, all values of 
USOUT should tend toward zero. The spikes in USOUT at 0.35, 0.40, 0.80, 
and 1.05 revolutions correlate with the impacts between roller and cage as 
shown in Figure 2. The first, second, and fourth events are the roller 
impacting the cage. This large impact force at the front of the roller 
tends to decelerate the roller in the path about the center of the bearing. 
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Figure 5. The Relative Slip Velocities Between the Rollers and the Inner Race Versus Time Measured in 
Revolutions of the Inner Race at 70,000 RPM 
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Figure 6. The Relative Slip Velocities Between the Rollers and the Outer Race Versus Time Measured in 
Revolutions of the Inner Race at 70,000 RPM 


thus giving the negative relative slip velocity. When the cage impacts 
the roller as in the third event, (0.80 revolutions) the impact force is at 
the rear of the roller. This tends to accelerate the roller in the path 
about the center of the bearing, thus giving a positive relative slip velocity, 

A measure of the overall skip in a high speed bearing is the deviation 
of the cage angular speed from the theoretical epicyclic cage angular speed. 
Figure 7 shows the variation of the cage angular spaed as a ratio of the 
theoretical epicyclic speed. This ratio hovers near 0.87 (a cage slip of 
13%) but still appears to be increasing, indicating that the new equilibrium 
has not yet been reached. The step changes correlate with the impact events 
described in the discussion of Figure 6, 

The trajectory of the cage center is a measurable quantity in a test 
of an experimental bearing. The dynamics of the cage are very important 
in the design of a high speed bearing as the cage is the component which 
provides the coupling effects among the rollers. Figure 8 shows the tra- 
jectory of the cage center in a coordinate frame whose dimensions are in 
units of the radial clearance between the cage and the inner race guide 
flange. 

Note that the cage has come in contact with the guide flange of the 
inner race at the bottom of the bearing near the dimensionless time of 
1,2 revolutions. This should tend to accelerate the cage, but the effect 
is not easily discernible from Figure 7, 

The cage-roller impact events presented in Figure 2 correlate with 
the abrupt changes in the cage center trajectory at the "knees" of the 
curve. The cage-roller impact exerts a radial force on the cage through 
the mechanism of rubbing traction between the cage and the roller in 
contact. 
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Figure 7. The Ratio of the Cage Angular Speed to the Theoretical Cage Epicyclic Angular Speed Versus Time 
Measured in Revolutions of the Inner Race at 70,000 RPM 
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Figure 8. The Trajectory of the Cage Center Normalized with Respect to the Radial 
Clearance Between the Cage and the Inner Race Guide Flange 
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Three Dimensional Simulation 


The three-dimensional simulation initial conditions and loading 
parameters represent the effects of a suddenly applied steady misalign- 
ment and a sudden switch to a steady shaft deflection in the downward 
vertical direction from the rotating shaft deflection condition. The 
maximum roller load stayed constant at 591 N. Figures 9 through 14 
indicate that the transient behavior due to the abrupt change in loading 
is completed by 0.20 revolutions. The general discussion of Figures 9 
through 14 and Figure 18 parallels the discussion for Figures 2 through 8. 

As it was for the two-dimensional simulation, the large values of FC 
(Fig, 9) correlate with abrupt changes in the USIN (Fig. 12), the USOUT 
(Fig. 13) and the cage velocity ratio (Fig. 14). The maximum impact force 
between a roller and the cage is (.171) (591) = 101 N (22.7 Ib^) . This 
force corresponds to a maximum Hertz stress betweeen the cage pocket and 

Q 

roller A of 3,26 x 10 Pa (41,500 psi) , 

The loading applied in the three-dimensional simulation is steady, 
both for the radial deflection of the inner race relative to the outer 
race and for the misalignment of the inner race relative to the outer race. 
The load zone does not move in this case, thus a roller requires a dimen- 
sionless time of 0,6 revolutions of the inner race to traverse it as shown 
in Figure 10. Note too that the amount of misalignment is very small and 
does not induce a load zone between the inner race and rollers opposite 
to the load zone induced by the shaft deflection. 

The magnitudes of the dimensionless traction forces, FIN and FOUT , 
in Figures 10 and 11 are comparable to those in Figures 3 and 4, The 
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Figure 9. The Dimensionless Normal Cage Forces Acting on the Rollers for the 3-D Simulation Versus Time 
Measured in Revolutions of the Inner Race at 70,000 RPM 
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Figure 10. The Dimensionless EHD Traction Forces on the Rollers at the Inner Race for the 3-D Simulation 
Versus Time Measured in Revolutions of the Inner Race at 70,000 RPM 



maximum value of the inner race/ roller traction forces is about 2.1 N 
(0.48 Ibj.) due to the low normal forces between the rollers and the 
inner race. The outer race/ roller traction forces reach values of (.05) 
(591) = 29.5 N (6.65 Ib^) , Figure 11. The plots of FOUT in Figure 11 
indicate that the dimensionless traction forces between the outer race 
and the roller hover about zero (0 to 0.01). The large variations in 
FOUT correlate to the rollers going through the load zone as shown in 
Figure 10. 

The relative slip velocities USIN and USOUT in Figures 12 and 13 
have the same characteristics as the corresponding plots in Figures 5 and 6. 
The values of USIN range from 12 to 16 m/s. The values of USOUT tend toward 
zero except for the spikes shown at 0.75, 0.80 and 1.0 revolutions. These 
events correlate with the cage to roller impacts shown in Figure 9. 

The plot of the cage velocity ratio versus dimensionless time in units 
of the revolutions of the inner race, Figure 14, remains close to 0.87. 

This plot is similar to Figure 7 in that both indicate a very small rate 
of increase in cage angular velocity and both show the sensitivity to cage- 
roller impacts. The sharp increase and decrease in the cage velocity ratio 
at 0.72 revolutions correlates with the cage/roller impacts shown in Fig- 
ure 9. 

The axial displacements of the rollers are plotted in Figure 15. The 
initial conditions for all these displacements were zero. The rollers 
which start to move axially are those which enter the load zone and are 
subjected to the roller skew moment induced by the misalignment. The 
rollers which pass through the total load zone, rollers B, A and G, are 
displaced the farthest away from the center, reaching three-tenths of the 
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Figure 11a. The Dimensionless EHD Traction Forces on the Rollers at the Outer Race for the 3-D Simulation 
Versus Time Measured in Revolutions of the Inner Race at 70,000 RPM 
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Figure 11b, The Dimensionless EHD Traction Force on Roller B at the Outer Race for the 3-D Simulation 
Versus Time Measured in Revolutions of the Inner Race at 70,000 RPM 
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Figure 12. The Relative Slip Velocities Between the Rollers and the Inner Race for the 3-D Simulation Versus 
Time Measured in Revolutions of the Inner Race at 70,000 RPM 
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Figure 13. The Relative Slip Velocities Between the Rollers and the Outer Race for the 3-D Simulation Versus 
Time Measured in Revolutions of the Inner Race at 70,000 RPM 
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Figure 14. The Ratio of the Cage Angular Speed to the Theoretical Cage Epicyclic Angular Speed for the 3-D 
Simulation Versus Time Measured in Revolutions of the Inner Race at 70,000 RPM 
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Figure 15. 


The Dimensionless Axial Displacement of the Rollers Versus Time Measured in Revolutions of the 
Inner Race at 70,000 RPM 




total axial clearance between the roller ends and the inner race 
guide flanges. 

The rollers are impacting the inner race guide flanges at dis- 
crete times during the 1.20 revolution period. As these events are of 
extremely short duration and are not continuous, it would be difficult 
to capture these data for plotting. A statement is printed by the com- 
puter program whenever the restoring torques due to the impact is non- 
zero, which permits the identification of the event time and magnitude. 

The skew angles of the rollers are plotted in Figure 16. The magni- 
tudes of the skew angles are small due to the very small axial clearance 
between the rollers and the inner race guide flanges. The skew angles 
fluctuate from positive to negative. This indicates that the rollers are 
exhibiting a "barrelling” motion as they move from side to side while 
orbiting in the bearing. 

The dimensionless skew moment plotted in Figure 17 is always positive 
for this loading as the effect of the misalignment only presents itself 
in the primary load zone. The skew moments correlate to the loaded rol- 
lers as inferred from Figure 10. The sharp irregular curves are a result 
of subtracting the skew moment effect of the smoothly varying inner race 
traction forces from the skew moment effect of the irregular outer race 
traction forces. Figure 11. 

The cage trajectory for the three dimensional simulation. Figure 18, 
has a different characteristic compared to the trajectory for the two di- 
mensional simulation. The trajectory in Figure 18 is essentially a straight 
line until the cage impacts the inner race guide flange and changes direc- 
tion near the bottom of the plot. The radial forces due to cage/roller 
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Figure 16. The Skew Angle of the Rollers Versus Time Measured in Revolutions of the Inner Race at 70,000 RPM 
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Figure 17a. The Dimensionless Skew Moment Acting on the Rollers Versus Time Measured in Revolutions of the 
Inner Race at 70,000 RPM 
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Figure 17b. The Dimensionless Skew Moment Acting on Roller B Versus Time Measured in Revolutions of the 
Inner Race at 70,000 RPM 
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Figure 18. The Trajectory of the Cage Center Normalized with Respect to the 
Radial Clearance Between the Cage and the Inner Race Guide Flange 
for the 3-D Simulation 


impact have a minor effect in this case as compared to Figure 8, 

This discussion has shown that misalignment leads to roller skewing, 
that the roller regularly impacts the guiding flange and that roller/ 
cage impacts can significantly affect the cage trajectory. In both the 
two-dimensional simulation and three-dimensional simulation, the slip 
velocities of the rollers at the inner race contact range between 12 and 
16 m/s. The slip velocities of the rollers at the outer race contact stay 
close to zero. These effects are the result of the rollers trying to achieve 
force and torque equilibrium subject to relatively large forces at the outer 
race and relatively small forces at the inner race. 

The results showing the effects of misalignment and roller/flange 
impacts are in agreement with Gupta [4], Although two different bearings 
and operating conditions are being compared between this study and [4] , 
the slip velocities at the inner and outer races appeared to be relatively 
constant whereas Gupta shows linearly varying slip velocities and a step 
change when load is removed from the contact. These trends observed in [4] 
do not agree with the trends observed in the present report. 

One future modification which could be made to the computer program is 
a more detailed model for roller/flange interaction which would consider 
flange layback angle. The capability of specifying lubricant temperature 
at several key locations in the bearing would enhance the usefulness of 
the program and make it compatible with "steady state" dynamic simulation 
programs such as CYBEAN [11]. 

Future studies with this computer program should include a transient 
dynamic analysis of a production bearing which would include a parametric 
study of such important parameters as cage material, lubricant type, roller/ 
cage clearances and diametral clearances. The effects of start up, speed 
changes, and shut down should also be included. 
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CONCLUDING REMAEKS 


A computer program for solving the nonlinear differential equations 
of motion of a high speed lightly loaded cylindrical roller bearing has 
been developed. This computer program has the capability of performing a 
two-dimensional or a three-dimensional simulation. Two example runs were 
presented to illustrate the use of the program and the plotted output. The 
computer program will be available through NASA. 
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APPENDIX A 


PROGRAM INPUT 

The input to the computer program may be classified into the cate- 
gories of roller bearing geometry, material constants, operating speeds, 
loading conditions, controls for the computations and output, and a set 
of initial conditions. The input data deck for a two dimensional program 
consists of 13+(NZ) cards, while the three dimensional program requires an 
input data deck of 13+2 (NZ) cards. The extra cards are needed to specify 
the initial conditions for the additional 3(NZ) variables. 

Card 1: NZ, NPTS, NTYPE 

Format: 315 

NZ is the number of rollers. 

NPTS is the number of slices in each roller which is used in the load 

distribution calculation. 

NTYPE is the code used to specify whether a two or three dimensional run 
will be made (2 or 3) . 

The program has been dimensioned to accommodate 30 rollers and a cage 
for the two dimensional problem or the three dimensional problem. If 
larger problems are to be run, the dimension statements will have to be 
changed . 

Card 2: El, E2, E3, Ul, U2, U3 

Format: 6F10.6 

El, E2, E3 are the moduli of elasticity of the inner race, roller and 
outer race respectively (psi) . 

Ul, U2, U3 are the Poisson’s ratios of the inner race, roller and outer 
race respectively. 
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Card 3: 


DRUACE, DRROLL, PD, ROLCLR, RADCOR, FLOOR 


Format; 

DRRACE 

DRROLL 

PD 

ROLCLR 

RADCOR 
FLOOR 
Card 4 : 
Format ; 

L 

RCWN 
ROLWGT 
ROLDEN 
Card 5 : 
Format : 
CGPOLR 
CGDIA 
CGWGT 
CGLGH 

CGCLR 

CGPOC 


6F10.6 

is the diameter of the inner race (in.), 
is the diameter of the roller (in.), 
is the diametral clearance (in.). 

is half the total clearance between the roller ends and the 

inner race guide flange (in.). 

is the corner radius of the roller (in.). 

is the corner radius of the inner race guide flange (in.). 

L, RCWN, ROLWGT, ROLDEN 
3F10.6, E20.14 

Is the axial length of the roller (in.), 
is the crown radius on the roller (in. ) . 
is the weight of a roller (Ib^) . 

2 4 

is the density of the roller material (Ib^ - sec /in ). 
CGPOLR, CGDIA, CGWGT, CGLGH, CGCLR, CGPOC, CGHT 
E20.14, 6F10.6 

2 

is the polar moment of inertia of the cage (Ib^- in - sec ) . 
is the bore diameter of the cage (in.), 
is the weight of the cage (Ib^) . 

is the face width of the cage next to one side of the inner 
race guide flange (in.). 

is the radial clearance between the cage and the inner race 
guide flange (in.). 

is the clearance between the roller and the cage pocket sides 
in the circumferential direction (in.). 
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CGHT 


is the radial thickness of the cage (in.) 


Card 6: 
Format : 
VISCCP 

OILDEN 
DECFUL 
Card 7 : 
Format : 
ANGDEF 

OMANDF 

SHDEF 

OMSHDF 

Card 8: 
Format : 
WMXKPM 
Card 9: 
Format : 

H 


VISCCP, OILDEN, DECFUL 
3F10.6 

is the absolute viscosity of the oil at operating temperature 
(centipoise) . 

is the density of the oil (lb /in ) . 

m 

is the ratio of the oil volume to the total volume in the bearing. 
ANGDEF, OMANDF, SHDEF, OMSHDF 
5F10.6 

is the total angular misalignment between the shaft and the 
outer race (radians) . 

is the control variable for steady angular deflection about the 
X axis (0), angular deflection rotating at shaft speed WS(+1), or 
steady angular deflection about the y axis (~1) . (See Fig. 1 [1]) 
is the total shaft linear deflection relative to the outer race 
support (in,), 

is the control variable for steady linear deflection in the +y 
direction (-1) or a linear deflection rotating at shaft speed 
WS (+1) . 

WMXRPM 

FlO.l 

is the maximum shaft speed (rev/min) . 

H, T, TEND, TSPACE 
4F10.8 

is the starting sten size for numerical integration (dimension- 
less time). Refer to eqn (73) Part I [1], Choose a value of 
0,0002, If it is too large or too small, the integration program, 
GEAR, will adjust the step size. 
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T 

TEND 

TSPACE 

Card 10: 
Format ; 
TPLT 

Card 11: 
Format : 
Y(l) 

Y(2) 

Y(3) 

Card 12: 
Format : 
Y(4) 

Y(5) 


is the initial value of dimensionless time 
is the final value of dimensionless time 

is the interval of dimensionless time at which the vector of 
dependent variables is printed. 

TPLT 

F10.8 

is the interval of dimensionless time at which the quantities 
of interest are placed into a file to be plotted. 

Y(l), Y(2), Y(3) 

3D20.10 

is the initial velocity of the cage in the x direction non- 
dimensionalized by the theoretical epicyclic speed (eqn. 3) and the 
radial clearance. 

is the initial displacement of the cage in the x direction in 
units of the radial clearance between the cage and the inner 
race guide flange. 

is the initial velocity of the cage in the y direction non- 
dimensionalized by the theoretical epicyclic speed and the 
radial clearance. 

Y(4), Y(5), Y(6) 

3D20.10 

is the initial displacement of the cage in the y direction in 
units of the radial clearance between the cage and the inner 
race guide flange. 

is the initial angular velocity of the cage nondimensionalized 
by the theoretical epicyclic speed. 
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Y(6) is the initial angular position of the cage (radian). 

Card 13 through Card 12+NZ: Y(J), Y(J+1), Y(J+2) 

Format: 3D20.10 

where J=7, 10, 13, ..., 4+3(NZ) 

Y(J) is the initial relative velocity between the i ’th roller 
[J=4+3i] and the cage nondimensionalized by the clearance 
(CGPOC) and the theoretical epicyclic speed. 

Y (J+1) is the initial relative displacement between the i ’ th roller 
and the cage nondimensionalized by the clearance (CGPOC) . 

Y(J+2) is the initial angular velocity of the roller about its own 
axis nondimensionalized by the theoretical epicyclic speed. 

Card 13+NZ; FORDIM, OMTHEO 
Format: 3D20.10 

FORDIM is the force F* used in determining the time scaling. It is 
set to zero and determined internally for new runs. When the 
program is restarted, this is the value utilized in the pre- 
vious run (Ib^) . 

OMTHEO is the theoretical angular velocity used in determining the 

time scaling. It is normally set to zero for new runs. When 
the program is restarted, this is the value utilized in the 
previous run (rad/sec). 

The data cards described above are for a two dimensional run (NTYPE=2) . 
If a three dimensional run were to be made, NZ data cards would have to 
be inserted between Card 12+NZ and Card 13+NZ given above. The other data 
cards would remain the same. 

Card 13+NZ through Card 12+2(NZ); Y(K), Y(K+1), Y(K+2) 
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Format: 3D20.10 


where K=3(NZ)+7, 3(NZ)+10, ...» 6(NZ)+4 

Y(K) is the initial angular velocity of skew of the i ’ th roller 

[K=A+3(NZ)+3i] . It is nondiraensionalized by the theoretical 
epicyclic speed. 

Y(K+1) is the initial angle of skew of the roller (radian). 

Y(K+2) is the initial axial displacement of the roller. It is non- 

dimensionalized on the total clearance between the roller ends 
and the guide flanges (2*R0LCLR). 

The card with FORDIM and OMTHEO is always the last card in the data 

deck. 
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APPENDIX B 


PROGRAM OUTPUT 

The printed output is written on unit 6 at an interval specified in 
the input data. The state variables and other data needed to restart the 
program are written onto a unit labeled JDATA. In this program JDATA is 
unit 10. Most of the interesting output data are not the dependent vari- 
ables themselves, but such quantities as the traction forces between the 
rollers and the races or the slip velocities between the rollers and the 
races. These data, which are to be plotted, are written on a unit labeled 
JPLT, and are also plotted by the program. In this program JPLT is unit 9. 

After the differential equation solver has numerically integrated the 
differential equations to TEND, the program will stop. The last data to be 
written on the unit JDATA are a set of card images containing all the in- 
put data on the first ten cards as described above with the exception of 
the ninth card where the entry for T is replaced by the final value of di- 
mensionless time attained by the program. The initial conditions as de- 
scribed above are replaced by the values of the dependent variables at the 
final value of dimensionless time. The last card image contains the values 
of FORDIM and OMTHEO used in the computation. 

To restart the program using these card images, a new value of TEND 
must be inserted such that TEND is greater than T. This new set of data 
allows the user to examine the output of the program in stages and, if it 
appears to be functioning properly, to allow the program to continue. 

The plot data stored on the unit labeled JPLT can then be trans- 
ferred to a program which controls the plotting routines. Every computer 
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has its own specific library of plotting- routines so the most general 
way of presenting the output data is to store it in a data file and adapt 
it to an available computer system and software. A subroutine, PLODAT, 
is used to plot the data on the University of Illinois CYBER 175. 

The data is stored in blocks on unit JPLT. The first line of each 
block has the counter value associated with that block, IPLT, the time 
in units of shaft revolutions, PTREV, the x and y positions of the center 
of the cage scaled on the radial clearance between the cage and the inner 
race guide flange, PXCG and PYCG, and the angular speed of the cage scaled 
on the theoretical cage epicyclic speed, PCGV. The format used for this 
line is (I 10, 7 E 10.4). 

The remaining data in the block is given, for each roller, as: 

The cage-roller interaction force on the roller scaled by F*. 

FC(J), J=1,NZ 
FORMAT (8E10.4) 

The traction force on the roller at the inner race contact scaled 
by F* . 


FIN(J), J=1,NZ 
FORMAT (8E10.4) 

The traction force on the roller at the outer race contact scaled 
by F*. 


FOUT(J), J=1,NZ 
FORMAT (8E10.4) 

The relative sliding velocity between the roller and the inner race 
in (m/sec). 
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USIN(J), J=1,NZ 
FORMAT (8E10.4) 

The relative sliding velocity between the roller and the inner race 
in (m/ sec) , 


USOUT(J), J=1,NZ 
FORMAT (8E10.4) 

The following data are also plotted for the three dimensional 
problem: 

The roller skew angle in units of radians, 

Y(J), J=3*NZ+8, 6*NZ+5, 3 
FORMAT (8E10.4) 

The roller axial displacement scaled on the total axial clearance 
between the roller and the cage pocket, 

Y(J), J=3*NX+9, 6*NZ+6, 3 
FORMAT (8E10.4) 

The skewing torque on the roller induced by non symmetrical EHD 

2 

traction forces on the roller scaled by the product I co , where I is the 

t e t 

transverse moment of inertia of the roller and m is the inverse of the 

e 

characteristic time associated with the scaling force F*. 

ROLMOM(J), J=1,NZ 
FORMAT (8E10.4) 
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